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Heat fluctuations over a time r in a non-equilibrium stationary state and in a transient 
state are studied for a simple system with deterministic and stochastic components: a Brow- 
nian particle dragged through a fluid by a harmonic potential which is moved with constant 
velocity. Using a Langevin equation, we find the exact Fourier transform of the distribution 
of these fluctuations for all r. By a saddle-point method we obtain analytical results for the 
inverse Fourier transform, which, for not too small r, agree very well with numerical results 
from a sampling method as well as from the fast Fourier transform algorithm. Due to the 
interaction of the deterministic part of the motion of the particle in the mechanical potential 
with the stochastic part of the motion caused by the fluid, the conventional heat fluctuation 
theorem is, for infinite and for finite r, replaced by an extended fluctuation theorem that 
differs noticeably and measurably from it. In particular, for large fluctuations, the ratio 
of the probability for absorption of heat (by the particle from the fluid) to the probability 
to supply heat (by the particle to the fluid) is much larger here than in the conventional 
fluctuation theorem. 

PACS numbers: 05.40.-a, 05.70.-a, 44.05. +e, 02.50.-r 



I. INTRODUCTION 

Knowledge of the behavior of heat fluctuations has an intrinsic value, especially for small systems 
such as nano-systems and bio- molecules, where the fluctuations are relatively large, but the wide- 
spread interest in fluctuation theorems (FTs) stems mostly from the fact that although there are 
few general results in non-equilibrium statistical mechanics, these theorems seem to provide some. 
The conventional FTs state that the probability distribution functionfPDF) ir for the average over 
a time r of a physical quantity A to have a value a, satisfies 0, 0,0,3,0,0 

7r((A) T = a;r) 

— — r « exp ar . 1) 

n{{A) T = -a;r) 

where (A) T denotes the time average of A and the dependence of the PDF on r has been explicitly 
indicated. The asymmetry in Eq. (^Q) is caused by some external field that is present in these models, 
such that a positive value of (A) T can correspond to a behavior "with the field" and a negative 
value to a behavior "against the field". Note that for large r, positive values for a become much 
more probable than negative ones. To be more concrete, in dynamical systems theory, the quantity 
A pertains to fluctuations in the "phase space contraction rate" Q, while in stochastic systems, 
A involves an "action functional" |6j . In specific cases, both have been connected with the heat or 
entropy production. Such a connection is based on the expression of the phase space contraction 
rate and the action functional, respectively, which typically have the form of a thermodynamic force 
times a current, divided by the temperature of the system^, i.e., the same form as the entropy 
production in Irreversible Thermodynamics. It is however not necessary to make the connection 
with the entropy production and in this paper, we will simply use the heat. 

There are in fact two different kinds of FTs: a transient (TFT) and a stationary state FT 
(SSFT). While the conventional SSFT only holds for sufficiently large (strictly infinite) times, the 
conventional TFT holds as an identity for all times [a]. Thus, for the SSFT, the ~ sign in Eq. (|T]) 
indicates the large r behavior, whereas for the TFT, it can be replaced by an equality sign|9|. 
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Apart from an (early) laboratory experiment on the SSFtQ, the FTs were restricted to theo- 
retical and simulation approaches: it was difficult to make laboratory experiments on macroscopic 
systems, since the large number of particles reduces all fluctuations enormously. Recently, Wang 



et al.\ll\ measured a TFT in the laboratory, by studying the motion of a single Brownian particle 



n part 

dragged through water by a laser-induced moving (confining) potential. But while in Ref. the 
entropy production (or heat) fluctuations over a time r were intended to be studied in a transient 
state, in fact the fluctuations in the work done on the system during that time were studied. These 
differ from the heat fluctuations due to the joint presence of the confining potential and the water, 
as we will explain more clearly later in the paper. While for this system the conventional SSFT 
and TFT do hold for the total work done on the system [l2T]. for heat fluctuations, very different 
FTs hold, in which the behavior of large heat deviations is notably and measurably different from 
the conventional FTs, as shown in a recent letter Here we present the full theory underlying 
the conclusions of Ref. 0] • 

The paper is organized as follows. Sec. [H] contains the basic theory. We present the model 
(Sec. Ill A|) . explain the difference between work and heat in the model pi Bj) . and treat the work- 
related SSFT and TFT fully plTty. For the heat fluctuations there is no expression for the PDF 
in terms of known functions, but its Fourier transform can be calculated exactly piD|) . Since this 
has no known exact inverse, we first calculate, in Sec. II I II the PDF of the fluctuations using two 
numerical methods: a sampling method as well as the algorithm of the fast Fourier transform to 
perform the inversion of the Fourier transform. The results of both methods agree well with each 
other and point to violations of the conventional SSFT and TFT. Next, this is substantiated by 
an analytic method developed Sec. IIVI based on the saddle-point method, which is shown by a 
comparison with the numerical methods, to work well when r is not too small. In Sec. ^ we show 
how this method can be used to obtain expressions for the extended heat FTs, both in the limit 
of r — ► oo, as well as for finite times. We conclude with a discussion of the results in Sec. I VII At 
the end, two appendices give some technical details used in the text. 

II. BASIC THEORY 
A. Model 

Theoretical descriptions of the experiment of Wang et al. 11] have been given even before the 
experiment by Mazonka and Jarzynski 14 ] , as well as recently by the authors [I2I H^ |. All were 



based on an overdamped Langevin equation for the Brownian particle, which reads: 

it = -r- 1 (x t -x t *) + a- 1 G. (2) 

Here Xf is the (three-dimensional) position of the Brownian particle at time t, a = 6ttt]R is the 
(Stokes) friction of the particle in water, rj is the shear viscosity of water and R is the radius of 
the particle. Furthermore, the relaxation time of the position of the particle in Eq. l(T|). r r . equals 
a/k, where k is the strength of the harmonic force. This force is derived from the potential 

£/( X ,t) = ^|x-x:| 2 , (3) 

where X j is the position of the minimum of the harmonic potential at time t, given in our case by 

x* = v*i, (4) 

(for t > 0) with v* the (constant) velocity of the motion of the potential through the water. 
Finally, the random force £t is taken to be Gaussian and delta function correlated in time: (Q) = 
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and (CtCs) = 2kBTa5(s — t)t. Here, the brackets denote an average over realizations, T is the 
temperature, ks is Boltzmann's constant and 1 denotes the 3x3 identity matrix. 

In the rest of this paper, we will use a reduced set of units. For our time unit we will use r r : 
thus wherever t appears one should read t/r r ; for our energy unit, we choose k^T and for our 
length unit, we use the width of the equilibrium PDF of Xt in the potential, i.e. \Jkj^Tjk. This 
change of units has been carried through consistently in the paper by simply setting a = 1, k = 1 
and ksT = 1. Then the Langevin equation above in Eq. (j2J takes the simple form 

± t = -(x t -v*i) + Ct, (5) 

where Ct satisfies 

(Ct) = 0, (6) 

{( t ( s ) = 25(t-s)t. (7) 

In Ref. ^3], the distribution of Xt and its time autocorrelation function were determined, under 
the assumption that the particle has been in the potential from time t = — oo up to t = and 
that during that time the potential did not move, so that the initial PDF at time t = is the 
equilibrium one. We found for the PDF of the position x of the Brownian particle at time t: 

p(x; t) = (2^)- 3 / 2 e-l x - v *^ 1+e " t )l 2 / 2 . (8) 

The resulting average position (xt) is therefore 

(x t > = v*(i - 1 + e-*). (9) 

Equations (JHJ) and © show that the system is not in equilibrium, but that, after a transient time 
period of order 1 (after which one may neglect the term e~ t ) the system becomes stationary in a 
co-moving frame — the PDF shifts to the right with a constant velocity v*. The mean position 
of the particle then becomes v*(t — 1), which is not equal to the position of the minimum of the 
potential, xjf = v*i. The difference — v* shows we are in a stationary state where the average 
harmonic force — ((x^) — v*t) = v* balances the friction force — v*. 

Below, we will need the correlations of the position at different times. For the time- 
autocorrelation function of x^, Ref. ^3] gives in current units: 

(x i2 x il )-(x i2 )(x tl )=e-l* 2 -' 1 ll. (10) 

B. Work versus heat fluctuations 

We first consider the fluctuations in the total work done on the system. This work is the total 
amount of energy put into the system. By the system, we mean here the particle in the potential 
U plus the water. At any time, the harmonic potential exerts a force Ft = — (xt — v*i) on the 
particle. Consequently, the particle exerts a reaction force —Ft on the potential. The potential 
has to be kept moving at a fixed velocity v*, for which an external force on it of magnitude Ft is 
required to balance the reaction force. Hence, the fluctuating work W T done on the system over a 
time interval from t to t + r is given by 

W T = [ + dt' v* • Fj/ = f + dt'w* ■ [-(x t > - vY)]. (11) 
Jt Jt 
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We next consider the heat Q T and its fluctuations, which in contrast to the total work W T , is 
only that part of the work which goes into the fluid, while the other part of the total work goes 
into the potential energy of the particle in the harmonic potential. Therefore 

Q T = W T - AU T , (12) 

where 

AU T = U(x t+T ,t + T)-U(x t ,t). (13) 

The potential energy U was given in Eq. ©15]. Note that in Eqs. I|ll |) -p3| ) we suppressed the t 
dependence. 

Before we discuss the fluctuations in the work and the heat, we want to say a few words on 
their average behavior. Using Eqs. © and (|lljl. we find that the average work in a time interval 
from t to t + t is given by 

(W T ) =w[t-(1- e- T )e~ l \ , (14) 

where 

w = \v*\ 2 (15) 

From Eq. (|14j). one sees that after a transient time t of 0(1), the second term within the square 
brackets can be neglected. Then, i.e., in the stationary state, the average work is wt, which means 
that the rate of work w done on the system is constant. To calculate the average heat, we use 
Eq. (|12|) and consider first the average potential energy 

(U(x t ,t)) = y"dxp(x,t)i|x-v*t| 2 = | + l w (l - e"*) 2 , (16) 

[by Eqs. © and (|T5J)]. from which we find, by Eq. (fT5]). that (AU T ) = 
w [(1 — e~ T )e~ t — \{1 — e~ 2r )e~ 2t ] . Combining this with Eqs. (fT2|) and (|T4"|) . one obtains 



( Qr )= w [ T + ^ 1 - e -2T )e -2ty (17) 

From this equation, we see that for the heat too, after a transient time of 0(1), the rate of heat 
production becomes a constant. Furthermore, this constant is the same as that for the work in 
Eq. (|14|) . So in the stationary state, on average all work is converted into heat. One might therefore 
be tempted to identify heat and work, at least in the stationary state. However, it will become 
clear that such identification is not possible for the fluctuations in work and heat. 



C. Distribution of work fluctuations 

Here we will briefly summarize the results for the work fluctuations from Ref. ^3] needed in this 
paper. The work in Eq. (Jll[) is defined as an integral over the path the particle. Since W T is linear 
in the path x^, and the x^ are Gaussian distributed, it follows that W T is Gaussian distributed as 
well. Hence, only the first two moments of the PDF need to be considered. 

For the work-related TFT, let us first consider the PDF for work done on the system during a 
transient state of time r , i.e., for an ensemble of transient trajectories of duration t, starting in 
equilibrium at t = 00, @]- The first moment can then be found from Eq. ()14|) . with t = 0, 



(W T ) T = w (t-1 + e- T ) , 



(18) 



5 



where subscript T denotes that the average is over transient trajectories. The second moment can 
be calculated as in Ref. using the relations in Eqs. (|8|)-(|11[). to be 

([W T -(W T }] 2 ) T = 2w(T-l + e- T ). (19) 

This is exactly twice the average in Eq. ()18|) . which turns out to be crucial. Given the mean 
m and the variance v of a general Gaussian distributed variable s, its PDF is given by -P(s) = 
exp[— (s — m) 2 /2v]/y/2irv, hence 

= e 2ms/v_ (2Q) 



If the variance is twice the mean, v = 2m, the right-hand side (r.h.s.) of this equation becomes e s , 
which is of a similar form as the FT in Eq. To apply this equation to the work fluctuations in 
a convenient way, we define p as the average rate of work in time r, W t /t, scaled by the average 
rate w of work in the stationary state: 

WT 

The transient PDF P^ V (W T ; r) of W T and the transient PDF tv^(p; r) of p are related by a constant 
Jacobian: 

<(;p;t) =wtP™{W t -t). (22) 

Superscript W denotes that these are PDFs related to work. Using Eq. (|2Uj) and the fact that for 
s = W T the variance is twice the mean, one finds the conventional TFT: 

= *MW T ] = exp[wrp]. (23) 

Note that the Jacobian drops out on the left hand side of this equation. Equation 1)23)1 coincides 
with Eq. flTJl. with {A) T = wp, and shows that the TFT holds for work fluctuations. In fact, it 
holds exactly ("as an identity" ja] ) for all times r. 

For the work-related SSFT, one should in principle look at a single (half-infinite) trajectory in 
the stationary state, and consider the statistics of work done in time intervals of length r along 
that trajectory. However, because of its stochastic nature, our system behaves ergodically and the 
distribution of the initial points of these intervals is simply the stationary one. Thus, the desired 
statistics can be determined by considering a single interval r of which the initial point is sampled 
from the stationary state. In addition, as shown in the previous sections, the system reaches a 
stationary state as t — > oo, in an exponential fashion. Hence we can generate the sampling of the 
initial point of the interval by considering the interval [t,t + r] for t — ► oo. In that limit, Eq. 1)141) 
gives 

(W T ) S = wr. (24) 

The subscript S denotes that the average is to be taken over trajectories in the stationary state, 
in the t — > oo limit just explained. The variance can be computed similarly, and the result is[l2T| 

([W T -(W T )] 2 ) s = 2w(T-l + e- T ). (25) 



Note that now the variance in Eq. (|25|) is not equal to twice the mean in Eq. 1)24)) . except when 
r — ► oo. We therefore expect that only for r — > oo, we have 

« exp[^ r ] = expkwp], (26) 

\ Pi T ) 
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where tt^ is the PDF of the time-averaged scaled work p in the stationary state. Since the r.h.s. 
goes to infinity in the r — > oo limit, Eq. (|26|) is not a well-defined result. A more precise form of 
the SSFT is: 



lim f^(p;r) = p, 



(27) 



where one defines the work fluctuation function as 



1 



In 



WT 



<(-p;r) 



(28) 



Indeed, using Eqs.J2 
fluctuations holds 12]. 



2D) and Eq. (j^ZJ follows, i.e., the conventional SSFT for the work 



D. Distribution of heat fluctuations 

The difficulty in considering the heat Q T instead of the work W T lies herein, that it is not linear 
in the position of the Brownian particle, through the contribution of U in Eqs. (|12|) - (jl3j) . which 
according to Eq. ® is quadratic in x 4 . As a result, the PDF of Q T , denoted by P® in the transient 
state and P® in the stationary state (with superscript Q for heat), will not be Gaussian and it 
does not suffice to calculate their first two moments. Nonetheless, the Fourier transforms of these 
PDFs, 

/oo 
dQ T e*rtTP?(Q T ;T), (29) 
-oo 

can be explicitly calculated, as will be shown below. Here, j stands either for T or for S, and this 
notation will be used throughout the paper. 

Consider first the transient case. We start by considering the joint PDF P^(W T , Axi, AX2; r) 
of W T , Axi and Ax2. Here, Axi and Ax2 are the positions of the particle relative to the position 
of the minimum of the potential U at times and t, respectively, i.e., 

Axi = x , (30) 
Ax 2 = x T - v*r. (31) 

Because W T , Axi and Ax 2 are all linear in the path xt, P^ is Gaussian, so that only its mean 
and variance are needed. Using Eqs. J5J), 1)12(1 and (fT3")l . the transient PDF P® of Q T is related to 
P^(WV,Axi,Ax 2 ;r) by 



P?{Qt\t) = JJJ dW r dAxidAx 2 P^(W T , Axi, Ax 2 ;r) 

x 5 (Q T - W T + ±[\Ax 2 \ 2 - |A Xl | 2 ]) . (32) 

The same can be done for the stationary case. The PDF P£ (W T , Axi, Ax 2 ; r) is defined as the 
infinite time limit of the joint PDF of W T , and the relative positions Axi and Ax 2 defined by 

Axi = xt - v*t (33) 
Ax 2 = xt +r -v*(t + r). (34) 

Like P£, Pg is Gaussian. From P s *, the stationary state PDF of Q T can be found from 

P s Q (Q r ;r) = JJJ d^ r dAxi(iAx 2 P*(W T ,Axi,Ax 2 ;r) 

x5(Q T -W T + l[\Ax 2 \ 2 -|A Xl | 2 ]) . (35) 
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The mean and variance of P* in Eqs. ()32|) and (|35j) are calculated in Add. lAl 

Equations (f32|) and (|35]> cannot be integrated explicitly, due to the quadratic nature of AU T , 
but their Fourier transforms can be obtained, as is seen when one combines Eq. ()32|) or ()35|) with 
Eq. (123) to yield 



(36) 



Because P* and e iq ^ WT d Ax2 l 2 l Ax il 2 )/ 2 l are both Gaussian, the integrals in Eq. (|36|) can be explicitly 
performed. The details of the calculation are given in App.[Bj with the result 



wq(i-q) 



where for j = T, A T = 1, while for j = 
the inverse Fourier transform 

P?{Qt\t) 



l + (l-e 



- T )] 



+ 



,-2 



r) g 2]3/2 



(37) 



S, A s = 0. The PDFs P? are related to the A(<?;t) by 



1 



dqe-^Pjiq-r). 



(38) 



Before we proceed with the evaluation of Pf(Q T ;r) from Eqs. (|3*7j) and (|38|). we note that if the 



conventional FT were to hold exactly, i.e., [cf. Eq. (JTJ with (A) n 



Qr/r] 



P?(Qt,t) 



P?(-Qt,t) 



(39) 



then, by the definition (|29|). we would have 

P j {q-T)=P j {i-q;T). (40) 

If Eq. (|4U|) does not hold, the fluctuation theorem in Eq. (|39|) cannot hold exactly for all r. Neither 
P T (g; r) nor P s (q; r) in Eq. (|37|) satisfy Eq. (|40j) . That means we can rule out the possibility of a 
TFT or a SSFT for heat fluctuations that holds as an identity for all r. 

The first term in the exponent in Eq. (|37|) does have the symmetry in Eq. (|40[) . and this term 
is the only one that grows with the time r. However, one cannot conclude from this that the FTs 
hold for large r, due to the singularities in Eq. (pTfj). In the rest of the paper, we will be concerned 
with whether Eq. (|39jl holds, both for transient and stationary states, for r — » oo, or, if it does not, 
what the deviations from this behavior are. To be precise, similar to the treatment in Sec. Ill CI we 
define the scaled heat fluctuations — denoted by the same symbol p as scaled work fluctuations — 
as 

P^— , (41) 

WT 

[cf. Eq. ((211) an d the remark below Eq. (JTZJ)]. Its PDF is 

(p; r ) = wrPf (wrp ; r) (42) 



WT 

2^ 



-iqwTp 



(43) 



[cf. Eq. (|22|)]. where Eq. (|3*5|) was used. Furthermore, we define the heat fluctuation function [cf. 
Eq. (25)] 



/, q Cp;t) = — In 



nf (p; r) 



-p;rj 



(44) 
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FIG. 1: The PDFs tt^ obtained by two numerical methods, for w = 2.0 and r = 1.0. Sampling method 
results (N s = 2, 000, 000, Sp = 0.2) are shown as crosses for tt® and as dots for 7rJ?, respectively. Fast Fourier 
Transform results (Sq = 4 ■ 10~ 4 , q max — 200) are shown as a thin dashed line for and as a bold dashed 
line for . The inset shows the corresponding fluctuation functions [Eq. I)44[l] and a solid straight line 
indicating the prediction of the conventional FT. 

and we investigate whether [cf. Eq. (|28|)] 

lim / Q (p;r) = p. (45) 

as is required for the conventional FT. 



III. NUMERICAL APPROACH 

Since no exact Fourier inverse of Eq. Q37[) is known, we first treat the problem of finding the 
PDF of heat fluctuations numerically. 

We used two numerical methods. The first is a sampling method which starts from the expres- 
sions in Eqs. and lf35|) that give the -P? in terms of the . From the Gaussian PDF P* with 
the mean and variance calculated in App.|Aj one draws many sets of values (W T , Axi, AX2), using 
a random number generator 0]. From each set, using Eqs. (fT^|) . ((T3|) and (|41[). p = (W t — AU t )/wt 
is calculated. From these p values we build a histogram by constructing bins of width 5p in an 
interval [—PmaxiPmax] and counting how many of the values fall into each bin. If the number of 
sample points N s is large and 5p is small, the histogram gives a good approximation for the PDF 
iTj. As a simple error estimate, one performs this procedure a few times with varying random 
seeds, and determines the spread in the results. 

The second numerical method is the standard Fast Fourier Transform algorithm [ill] applied to 
the inverse Fourier transform. This algorithm takes a discrete set of values of the function Pj(q) 
for q in an interval [—q ma x, Qmax] with equal spacing 5q between the values, and returns values of 
the (inverse) transform, i.e., of the function Pj 5 , on a reciprocal grid of values for Q T . If 5q is small 
enough, and q max is large enough, this can be used to obtain a good approximation for the (inverse) 
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1 2 3 1 2 3 

P P 



FIG. 2: Numerical results for the fluctuation functions fj for the stationary (j = S) and the transient case 
(j = T), for (a) w = l,r = 4, (b) w = 2,r = 2.5, (c) to = 3,r = 1.3333 and (d) w = 4,r = 1. Sampling 
method results (iV s = 1.2 • 10 9 , Sp — 0.2) are shown as thin crosses for and as bold dots for f§. Fast 
Fourier transform results (Sq = 4 ■ 10 -4 , <7 ma :r = 200) are shown as the dashed thin curves for and as 
dashed bold curves for f$ . The solid straight line indicates the prediction of the conventional FT. 

Fourier transform. An error estimate can be found by varying 5q and q m ax an d observing to what 
extent the results have converged. In all plots in this paper, these errors would be unobservable 
and are therefore not plotted. Note that once P® is known, by Eq. ()42|) we can obtain 7r? as well. 

The results of these two methods for the PDFs ir® and for the parameter values w = 2.0 and 
r = 1.0 have been plotted in Fig. ^ as a function of p, and the corresponding fluctuation functions 
/x and f® are plotted in the inset in that figure. The two numerical methods agree very well, but 
do not agree at all with the straight line with slope 1 of the conventional FT [i.e., p cf. Eq. (|45[)] . 
which is also drawn in the inset in Fig. ^ neither in the stationary nor in the transient case. 

To investigate the discrepancies with the conventional FTs further, we have explored a range 
of values for the parameters w and t, using the two numerical methods. Some results are shown 
in Fig. [21 One sees that for large values for p, deviations of f®(p; r) from the straight line p of the 
conventional FT are generic. 

Furthermore, Fig. |21 shows (for w = 2.0, but the same holds for other w values) that for large r, 
the straight line is approached for small p, both in the transient and the stationary case, although 
the approach is considerably slower for the former. 

IV. ANALYTIC APPROACH SADDLE POINT METHOD 

The numerical methods are not suitable to give reliable results for the tails of the PDFs for 
large r. To get around this, we developed an asymptotic analytic approach. 

The starting point of the approach is to notice that the exponent in Eq. (|H7|) grows linearly 
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FIG. 3: Numerical results for the r dependence of the fluctuation functions for w = 2.0 for (a) the 
transient (j = T) and (b) the stationary case (j = S). Sampling results (N s = 1.2 ■ 10 9 , Sp = 0.2) are shown 
as symbols: (*) r = 0.5, (o) r = 1, (•) r = 2, (□) t = 4. Fast Fourier Transform results (<5<7 = 4 • 10~ 4 , 
<Zmaa = 200) are shown as the dashed curves (the r values of these curves are the same as those of the 
coinciding sampling points). The solid straight line is the conventional FT result. 



with t. Thus, the 7r?, given by Eq. (|43[). can be written as 

tt?(p;t) = — / dqe- WT ^ q ri +i « p \ (46) 

J 27r y.^ 



where 



e i(</^) =-^M^(9;r)] (47) 



is of order r° [see Eq. (1371)1. The presence of a large parameter r in Eq. (|46[) makes it suitable for 
the saddle-point method |l7l|. which we will briefly describe now. 

To get a good approximation of the integral in Eq. ()46[) . one considers 

h(q) = -w[e j (q;T)+iqp] (48) 

as a function of a complex- valued q (where the dependences of h on j, p and r are suppressed). 
Then, one determines its saddle point, i.e., the complex number q = q* for which dh/dq = 0. 
Through this point q* in the complex plane lies a path of steepest descent S along which Im h is 
constant and Re h attains a maximum at q* (l^ |. Next, one continuously deforms the original path 
of integration R (here the real axis) to this path of steepest descent S without crossing singularities. 
The integral over S is for large r dominated by the (small) segment on S around the saddle point 
q* and can be expanded in inverse powers of r by Taylor expanding the function h on S around 
the saddle point. For the leading term one uses a second order Taylor expansion and finds [l7T| 



3 / 2-7T 

dqe rh ® = J^JF^ e ^*)+-[l + 0^)]. (49) 
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Here, a is the angle between the direction of the path S when it traverses the saddle point, and 
the real axis. Furthermore, h"(q*) is the second derivative of h with respect to q at the saddle 
point q*, which is assumed to be non-zero. The correction term in Eq. Q49|) is 0(r _1 ), though that 
strictly applies only when the function h does not depend on r. In our case it does, and we will 
see that consequently the correction terms can then be 0(t -1 / 2 ), although the leading behavior in 
Eq. JUJ) is not affected. 

The form of h(q) in Eq. (|48|) has the following simplifying consequences for the expression in 
Eq. (|49j). Firstly, the equation for the saddle point, dh/dq\ q * = 0, becomes 



e 



>(q*- T ) = -ip, (50) 



where e'- is the derivative of ej with respect to q. Secondly, h"(q*) in Eq. (|4*9")) is equal to —we"(q*). 
This second derivative can be further simplified by noticing that Eq. (|50|) holds for all p, so we can 
take the derivative with respect to p on both sides to find e'-{q*; r)(dq* / 'dp) = —i. Hence 



\h"(q*)\ = w\e"(q*;T)\ = w 



dq* 



dp 



(51) 



Thirdly, the approximation is for a PDF, which should be real and positive. Thus, if the saddle- 
point method is to be a consistent calculation of a PDF, we must have 

a = 0. (52) 

We will prove later that this is indeed the case. Using Eqs. (|48|) . (|51|) and (|52|) . we apply the saddle 
point approximation Eq. (|49|) to the integral in Eq. ()46|) to find 



7if (p;r) ~ 



/ WT 


dq* 


J 2^ 


dp 



e -WT[e j {q*;T)+ipq*\ _ ^) 



Here, the symbol ~ is for convenience used to denote the asymptotic behavior, instead of explicitly 
denoting the correction terms as 0{t~ 1 ) as we did in Eq. (|49|) . 

As mentioned above Eq. (J49|). this method gives an asymptotic expansion in inverse powers 
of t. Now if r is large enough so that terms of relative order r _1 can be neglected, then terms 
which are exponentially small in r can also be neglected jli|. Therefore, we use in the definition of 
the exponents ej in Eq. (|47l) the expression of the Fourier transform Pj(q;r) in Eq. (|37j) without 
exponential terms, so that[19j 

t \ c \ V + 3A jg 31n(l + g 2 ) 
ej ( g] r) = -g(, -q)- ^ - - + 2wr (54) 

When substituted in Eq. (|50j). we obtain a fourth order polynomial equation for the saddle points, 
so that there are four saddle points. To know which of the four saddle points to use, one needs to 
determine their location, the paths of steepest descent that traverse them, and to see if our initial 
integration line (the real axis) can be deformed to lie on one of these paths without crossing any 
singularities. 

For three typical values for p we have depicted the four saddle points of h for the stationary 
state in Fig. |1J calculated analytically by Mathematical^ which uses Ferrari's method for solving 
the quartic equation We see in Fig. Hj that for all p there is one saddle point (A) on the 
imaginary axis above i, one in between —i and i (B), and two (C and D) who are either both 
purely imaginary but below —i, or are complex with the same imaginary part (also below —i), and 
with opposite real parts. In Fig. ^Jthere are a pole at q = —i, and branch cuts which we have taken 
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FIG. 4: Geometry in the complex g-plane of the saddle points, singularities and paths of steepest descent 
of h(q) [Eqs. iji^ and for the stationary state [Aj = 0], for w = 1, r = 4, and p = —2, 0, and 4.25, 
respectively. In all graphs, solid dots depict the position of the saddle points A, B, C and D. Wiggly curves 
give the branch cuts ending in the branch points at q = ±i which are indicated with a thin horizontal line. 
The branch point at —i is also a pole. The dashed line R gives the original line of integration, i.e., the real 
axis. Solid lines labeled Sa, Sb, Sq and Sb are paths of steepest descent though the saddle point, on which 
Re h attains a maximum. 

along the imaginary axis, one from i upward and one from —i downward. These cuts come from 
the logarithmic term in Eq. (|54|) . Finally and most importantly, Fig. |1] shows the paths of steepest 
descent Sa, Sb, <S*c and Sd through the saddle points A, B, C, and D. These paths were found by 
having Mathematica solve for lines of constant Im/i, where the values of the constant for Imh are 
chosen to match the values of Im h at the four saddle points. 

Knowing now the geometry of h in the complex plane, we can determine which saddle point to 
use. One sees that for all p the path Sb through saddle point B is a possible deformation of the 
real axis: each point on the real axis is simply shifted up or down until it lies on this path and 
no poles or other singularities are crossed in doing so. Strictly speaking, the deformation changes 
the end points because they get shifted by an imaginary amount ia [where a can be shown to be 
(1 — p)/2], but one can add portions ia to the integration path from — oo to — oo + ia and from 
oo + ia to oo, and as the integrand vanishes on these portions, the integral taken over R and over 
Sb have the same value. This is not true for Sa, Sq and Sb- 

The graphs in Fig. suggest that the path Sb goes through B horizontally, i.e., that a = 
as Eq. (|5*2)) conjectured. To prove this, notice that for purely imaginary values for q, i.e., q = iX, 
h{i\) = — w[e(iX]r) — Xp] is real if —1 < A < 1 [cf. Eq. (|54|)]. Hence, the imaginary axis from 
—i to i is a line of constant imaginary part of h. In the saddle point, lines of constant imaginary 
part cross perpendicularly[12|. Thus in saddle point B, the curve of steepest descent is seen to be 
perpendicular to the imaginary axis, i.e., in the direction of the real axis, so that indeed a = 0. 

We remark that the above analysis for the stationary case can be carried through also for the 
transient case, with only a slight change in the saddle points. 

Our saddle-point method now consists of the following. We use Mathematica to determine 
analytically the purely imaginary solution q* of Eqs. (|5()|) using (|54j) with \q*\ < 1, as a function 
of p. Given that solution, the PDF 7r? is explicitly calculated using Eqs. (|53|) and (J54[). Once 7r? 
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FIG. 5: Comparison of the saddle-point approximation with the numerically obtained solution, for w = 0.5 
and (a) r = 2.5 and (b) r = 5. The solid thin line is the transient PDF tt®(p;t) as calculated using the 
saddle-point approximation of Sec. IIVI and the dashed thin curve is the same PDF calculated using the 
numerical inverse Fourier transform of Sec. Hill (So = 2 • 1(T 4 , q max = 2 18 5q = 52.4288). The solid bold line 
is the stationary PDF tt® (p; r) calculated using the saddle-point approximation, and the bold dashed curve 
is that same PDF obtained using the numerical inverse Fourier transform. 

is known, we determine the fluctuation functions from Eq. (|44|) . We can then investigate the 
FTs, even for finite r, something which is usually not possible, but is doable in this case because 
the Fourier transform is known exactly for all time. 

To see how good the saddle-point expression for 7r? in Eq. (|53|) is, we compare it with the 
numerical results from Sec. IIII1 specifically, the numerically determined Fourier inverse. This is 
shown in Fig. both for the transient and for the stationary case. In Fig. Ufa), we chose r = 2.5 
and w = 0.5, i.e., not too large, and we are therefore still able to see the discrepancies between the 
saddle-point result and the inverse Fourier transform for 7r?, which are about 15% near the peak, 
and just a few percents in the tails. In Fig. |SJb), we chose r = 5.0, which is still not too large, but 
we now see that the saddle-point method has become quite accurate, up to about 7% in the peak 
of the distribution, while in the tails one can hardly distinguish the results from the two methods. 
As r gets larger, the results from these two methods approach each other even more and plotting 
them together as in Fig. would show basically two curves on top of the other. 

While the numerical methods in Sec. IIIII could give us the large r behavior of f®(p;T) only 
for small p values, using the saddle-point approximation, we can now determine the behavior of 
ff{P'j T ) with increasing r for the full range of p values. The results are plotted in Fig. and show 
that both for the transient and for the stationary case, the fluctuation function approaches 
the conventional FT only for small p values. In contrast, for large p values (p > 3), a completely 
different limit emerges as r — > oo, one where ff{p, T ) appears to approach 2, both for the transient 
and for the stationary case. The exact form of this r — > oo limit of ffip, t) is given by Eq. I|68j) in 
Sec. El 

V. ANALYTIC EXPRESSIONS FOR EXTENDED FLUCTUATION THEOREMS 

Because we do not find it satisfactory to have the analytical result of the saddle-point method 
only in the memory of Mathematical, we will determine in this section explicitly the leading asymp- 
totic behavior of 7i\? and corrections as r gets large, using Eqs. (|50j) . ([53]) and (|54|) . 
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FIG. 6: Results for the large r behavior of (a) f®(p; r) and (b) f®(p] T ) f° r w = 1, as found from the saddle 
point approximation. The figures show four lines corresponding to result for r = 3, r = 10, r = 100 and 
t = 1000, as well as a solid black straight line (p) that gives the behavior expected from the conventional 
FT. 



To calculate the saddle-point approximation for ir® given in Eq. ((55]) . we need to find the 
saddle-point q* first. With ej from Eq. (|54j) . Eq. (j5()[) for q* takes the form 



4<f 2 (3-2ig*) + 3Aj 3iq* 

2r(i + q*) 2 u;r(l+^* 2 ' 



1 + + ^-7—3^ 1 " — T, ^ = P- (55) 



This equation always has four (complex) solutions as explained below Eq. (|54j) . Among these four 
solutions, one is always purely imaginary and lies between —i and i, which is the one we need, as 
it corresponds to saddle-point B (see previous section). In order to find this solution, we try the 
expansion 

q * =q * + Sl + (r- 2 ), (56) 
r 

where q$ is given by the solution of Eq. (|55|) with r — > oo, i.e., 1 + 2iq* = p. Hence 

q* =i(l-p)/2. (57) 

Since q* lies between — i and i, the above zeroth order solution is only appropriate if it lies in that 
interval, i.e., if —1 < p < 3. We will refer to this as case (a), which will be dealt with first. The 
cases (b) p < — 1 and (c) p > 3 will be treated later. 

(a) For — 1 < p < 3, Eq. ()57j) . when substituted into Eq. (J53f) using ej from Eq. yields a 
zeroth order solution for 7r?(p; r). Because (?q is the maximum of h in Eq. (|49|) . there is no need to 
calculate q\, as far as the exponent h is concerned |22j. Furthermore, in the prefactor in Eq. (|53|) . 
dq^/dp is of order 1, so adding r~ l dq\/dp gives a correction of relative order 0(r _1 ), i.e. of the 
same order as the correction terms in Eq. (|49|) . So there is no point to work out q* to higher orders 
than q$. Substituting the expression for g^, given in Eq. I|57|) . for q* in Eq. (|53|) and using Eq. (|54[) . 
yields 

[ ^— „[ r(1 _ p) . + ffi ; li { ( 1 _ p)a _s Aj} ] / 4 

*? ( ^> ~ [(3 +,)]■/» for-l<P<3, (58) 

with correction terms that are relatively 0(r _1 ). 
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Outside of the range — 1 < p < 3, the zeroth order Ansatz in Eq. (|5T|) fails to produce a solution 
for q* between —i and +i. The explanation is that the expression for q^ in Eq. (|57[) was based 
on neglecting the last two terms on the left-hand side (l.h.s) of Eq. (|55jl. which clearly does not 
work for p > 3 or p < —1. However, the only way that these terms can become comparable to the 
first term in that equation, is if their denominators become 0(1) instead of O(t) as they appear 
to be. For this another Ansatz than Eq. (|56[) is necessary. As there are two different denominators 
in Eq. (|55[) . there are two cases. In the first case, the denominator of the last term on the l.h.s. of 
Eq. (jSSJ) is taken to be 0(1), i.e., r(l + q* 2 ) = 0(1), or q* = ±i + 0(r _1 ). The sign ambiguity is 
lifted by noting that for q* = — i + 0(t _1 ), the second term in Eq. (|55|) would diverge oc r, so only 



q* = +i + 0{t~ 1 ) (59) 

is correct. In the second case, the denominator of the second term on the l.h.s. of Eq. (|55|) is taken 
to be 0(1), i.e. r(i + q*) 2 = O(l), or 

q * = - i + 0{T- 1 / 2 ). (60) 
(b) We first consider the case q* = i + 0(r _1 ) of Eq. (|59l) . which we write for convenience asj2^| 



i + — + 0{t- 2 ). (61) 



WT 

Substituting Eq. (|61jl into Eq. (|55|) . and keeping only 0(1) terms, one gets — 1 + j- = p, so 

« = 2fcW aIld 

3?' 

q*=i + f- - + 0(r- 2 ). (62) 

y 2wr(p+l) v ; v ; 

Therefore, for p < —1, this is an approximation of the solution q* below i. Using Eqs. (|53|) . ()54l) 
and (l62l). one finds 



7r?(p;r) ~ J^JLM±A e -^[-r P +i-|A,] + | for p < -1. (63) 

J V 36-7T 

Again, correction terms to this are of relative O^^ 1 ). 

(c) Finally, we consider the case q* = — i + 0(t -1 / 2 ) of Eq. (|6()j) and write q* = — i + ^= + ^. 
Note that in principle, we needed to include two correction terms here to get a similar level of 
approximation as in Eqs. (|56|) and (|61|), However, as the sum of the first two terms (i.e., —i + 
«&/ \/r) is to leading order the maximum of the exponent, we need not calculate cj^. Substituting 
q* = — i + ibT~ l l 2 into Eq. (|55[). gathering the 0(1) terms and using that = Aj, we find 
b = ±(2 — Aj)/y/2(p — 3). We need the positive solution, so that is above — i, i.e., 

q * = -i + i^^ + 0(T-% (64) 

Note that we need p > 3 for this to be purely imaginary, as it should be. Combining Eqs. (|53[). 
(J5H) and (jMJ) gives 



WT 2 



(p-2)r-(2-A i ) v /2(p-3)r+ 



40-12p+3A -(p-4) 



"?<" T ' ~ 1/ 16 *(2-A, V e 1 for p> 3, (05) 
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where A 2 = Aj has been used to simplify the expression. Due to the different form of the expansion 
of q* here compared to the cases (a) and (b), this result is valid up to corrections of relative order 
0{t- 1 / 2 ). 

Note that we now have the asymptotic forms for 7r? in the three separate regions, p < —1, 
— 1 < p < 3 and p > 3 in the Eqs. ()63j) . ()58j) and ()65|) respectively. This shows that the behavior 
in the center of the PDF is Gaussian-like, while in the tails it is exponential. 

To investigate the FTs, we will now consider the fluctuation functions For that, we have to 
take the above asymptotic expression for large but finite r to calculate of Eq. (|44jl (after which 



one can take the limit r — > oo). In this way, one finds for the large but finite r behavior of /• :[2 



/>;t) 



p + 



6A> 



3 In 

3— p — 2p 



r(9 -p 2 ) 



2wt 



+ 0(t- 



for < p < 1 



P 



(l-p) 2 /4- — + ^ + 0(r- 2 ) 
tor 



r,(p) 



(66) 



2 + (2 - Aj 



2(p - 3) 



for 1 < p < 3 
InT 



-1> 



2wt 



+ 0(t 
for p > 3, 



where 



In [w 2 (3 - p) 3 (l + pf(p - l)/576] 



3 w(p+l)(3Aj 
2 + 4(3 



2p 2 + 8p - 10) 



p) 



(67) 



To obtain Eq. we have used Eq. for < p < 1, and we have combined Eqs. (jSHj) and ljfi3*|) 
for 1 < p < 3, and Eqs. (|63|) and (|65[1 for p > 3. Only the case p > is given by Eq. ()66|) . For 
p < 0, one can use that by the definition Eq. (|44jl. /?(p; r) is an odd function of p. 

We remark that as we approach the limits of the range of p values for which each of the 
expressions for the PDF 7i\?, [i.e. Eqs. (|63|) . and (|65j) ] holds, their r.h.s.'s diverge. Similarly, 
the r.h.s.'s of Eq. (|66|) diverge there. These divergences are an artifact of the expansion in inverse 
powers of t, as can be seen because such divergences are not observed in Sec. IIVI (i.e., in Figs. El 
andlH. In fact, for any r, Mathematica shows us that q* is a continuous and differentiable function 
of p|2l| (with a discontinuity in the derivative dq* /dp only as r — > oo), leading to a continuously 
varying 7r?(p; r) with p. Furthermore, considering the regions of p values in which Eq. (|66|) differs 
(due to these divergences) noticeably from the full saddle-point approximation, one finds from 
Mathematica that these regions shrink to zero as r 

For r — > oo, the results in Eq. (|66|) become 



oo. 



lim / Q (p;r) 

t — >no J 



p 
p 
2 



(p-l) 2 /4 



for < p < 1 
for 1 < p < 3 
for p > 3 



(68) 



This is the extension of the conventional FT as it holds for infinite time. Comparing with Eq. (|45|). 
it is clear that the conventional FT only holds for not too large fluctuations, i.e. |p| < 1. For larger 
p, there is first a quadratic deviation after which the ratio between the probability for positive 
and for negative fluctuations becomes a constant. In contrast, for the conventional FT, f9 keeps 
increasing linearly with the size of the fluctuation p. 

Equation (|66[) gives the extension of the fluctuation theorem for finite r. The deviations from 
the infinite r results are of a different character for small (|p| < 3) and large (|p| > 3) fluctuations. 
The ff(p;r) for large positive fluctuations (p > 3) approaches 2 with corrections that scale as 
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l/yjr. But whereas in the infinite r limit in Eq. (|68|). the function no longer grows above 2 beyond 
p = 3, in the Eq. (|66|) for finite r this function keeps growing to leading order as a square root, i.e., 
f?(P'i T ) ~ V(P ~ 3)/ T - The prefactor of (p — 3)/r depends on whether we consider the transient 
or the stationary case — it is multiplied by \[2 for the former and by 2y/2 for the latter — but it 
does not depend on w. For small fluctuations, on the other hand, which might be more relevant 
for many applications, we see corrections to the infinite r limit of 0(r _1 ), which do depend on w. 
The value of w determines whether the curve starts around p = above or below the line with 
slope 1, and thus also whether the slope of 1 is approached from below or from above. In fact, 
expanding Eq. (|66j) around p = 0, we get 



2 /4-3A,- 1 
1 + - - 

r V 9 w 



p + 0(p 2 



(69) 



This shows that the critical value for w is w c = 9/(4 — 3Aj), i.e., w c = 9 for the transient and 
w c = 9/4 for the stationary case. For w = w c , the slope (near p = 0) is 1 for finite times (up to 
corrections of order r -2 ), while the slope of 1 is approached for large r from above for w > w c , 
and from below for w < w c . In contrast, for the conventional TFT, the slope is 1 for all r, while 
for the conventional SSFT, at least for the work done on the system, the slope always approaches 
1 from above irrespective of w with increasing r. 



VI. DISCUSSION 

1. This paper treats fluctuations in the heat developed in a system of a Brownian particle 
in water, confined by a harmonic potential, which moves at constant velocity through the fluid, 
dragging the Brownian particle with it. The theory of heat fluctuations developed in this system was 
based on an overdamped Langevin equation for the position of the particle. This theory required 
a far more sophisticated analysis than was used in the previous paper[lij for work fluctuations. It 
should be mentioned that some of this same sophistication is also found in the work of Farago for 
a quantity different from both the work and the heat (2^. 

2. In essence, our theory deals with the fluctuations of the quantities occurring in the first law 
of thermodynamics, i.e., work, heat and internal energy. The energy balance for the system is 

Q T = W T — AU T . (70) 

Here W T is the total work done on the system during a time r, Q T is the heat produced by the 
Brownian particle in water in that time, and AU T is the difference in potential energy of the particle 
in the harmonic potential in the same time interval. Equation (|7()|) . which is basically the first law, 
can be applied both to averages as well as to fluctuations, because it expresses energy conservation, 
which holds both macroscopically and microscopically. 

3. The theory gives extensions of the conventional SSFT and TFT in Eq. (|66|) . In the limit 
t — > 00, leading to Eq. (|68|) . this new theorem coincides with the conventional TFT and SSFT 
only for small fluctuations p < 1 (as Rey-Bellet and Thomas also found for a different system[26]), 
while for larger fluctuations the behavior is completely different from the conventional ones. We 
will now explain the qualitative behavior of work and heat fluctuations. For that, it is useful first 
to consider, in point 4 below, the system in equilibrium, i.e., in the situation in which the harmonic 
confining potential does not move. In point 5, we then discuss the non-equilibrium case. 

4. In equilibrium, work, potential energy and heat behave as follows, (i) Work: There is no 
displacement and hence no work is done. Therefore the PDF of W T is a delta-function: P(W T ) = 
5(W T ). (ii) Potential: The PDF of the potential energy of the particle in the harmonic potential, 
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U, is given by a Boltzmann factor ~ exp[— U] {ksT = 1), since the particle in the harmonic 
potential can be seen as a subsystem of a larger one. Prom this, one sees that the PDF of AU T 
behaves for large r similarly in its tails, i.e., P(AU T ) ~ exp(— | AU T \). (iii) Heat: For the PDF 
of the heat fluctuation, using Eq. (|70|) with W T = 0, we see that Q T also has exponential tails: 
P(Q T ) ~ exp(-|Q T |). 

5. We now turn to the non-equilibrium, stationary state case, (i) Work: the PDF of W T is now 
a Oaussian[l3]. (ii) Potential: The PDF of AU T is expected to have still exponential tails, at least 
near equilibrium, i.e., P(AU T ) ~ exp(— | AU T \). (iii) Heat: The effect of the interplay between W T 
and AU T on the behavior of heat fluctuations can be deduced using Eq. (j70|) . We have to separately 
consider small and large fluctuations. For large r, W T and Q T grow on average linearly in time, 
while AU T stays of O(l). Hence, for small fluctuations of these quantities (near their averages) one 
can neglect AU T in Eq. (J7UJ), so that Q T ~ W T , and the behavior of Q T is very similar behavior 
to that of W T , i.e. Gaussian-like. On the other hand, when a large fluctuation of Q T occurs, it is 
less likely to be due to a large fluctuation of W T than to a large fluctuation of AU T . This is so 
because the tails of the Gaussian PDF for W T are much smaller than the exponential tails of the 
PDF for AU T . As a result, W T will be near its average while AU T is large. For the sake of the 
argument, we put W T equal to its average, {W T ), which coincides with (Q T ) for large r. Hence by 
Eq. (231), Qt ~ (W T ) - AU T = (Q T ) - AU T . Using that the PDF for AU T behaves as exp(-|A[/ T |) 
in its tails, we see that the PDF for Q T has tails of the form exp(— \Q T — (Q T )\), in agreement with 
Eqs. (Pjl and lfH5j) . 

6. We saw in Ref. ^3] that the work obeys the conventional FT in the limit r — > oo. Heat and 
work fluctuations are expected to behave similarly for small fluctuations (cf. point 5), and therefore 
the conventional FT is obeyed also by the heat fluctuations for r — > oo for small enough fluctuations. 
However, for large fluctuations, we get a different behavior. Using P(Q T ) ~ exp(— \Q T — (Q T )|), the 
fluctuation function /? which is defined by Eq. (|44jl and can be written as ln[P(Q T ) / P(—Q T )]/ '(Q T ) 
becomes (— \Q T — (Q T ) | + \Q T + (Q T ) \)/{Q T ) = 2. This explains qualitatively the behavior expressed 
in Eq. fify . 

7. The symmetry relation in Eq. (|4()|1 is very reminiscent of the one used by Lebowitz and 
SpohnjS] in their work on the fluctuation theorem (i.e., e(A) = e(l — A)), although their method 
of large deviation theory allows only a treatment of the behavior of P® and /? for r — ► oo. The 
precise connection between their models (and the conventional FT they find), and our model (and 
the extended FTs) is not clear. 

8. We derived the extended infinite-r FT of Eq. (|68j) already in Ref. 13] using large deviation 
theory. In fact, the saddle-point method reduces to the large deviation theory of that paper in the 
limit t — > oo. We will not prove this here, but remark that if one takes r — > oo in Eq. (|54[) . then 
one gets —q{i — q). Setting q = iX, this becomes A(l — A), which is the form of the quantity e(A) 
used in the large deviation theory in Ref. for |A| < 1. The expression breaks down when the 
correction terms to —q(i — q) in Eq. (|54jl become infinite, i.e. at A = ±1. This restriction of |A| < 1 
was also crucial in obtaining the extended FT in Ref. [13j. We remark that likewise, in Fig. 01 the 
singularities at q = ±i restrict the saddle point to the region \q\ < 1 as well , and that this in turn 
leads to exponential tails of P(Q T ), which finally give rise to the extensions of the heat FT. 

9. One of the important and striking results of the extended FTs is that the probability ratio 
for negative to positive fluctuations in the heat production by the Brownian particle is much larger 
than that given by the conventional FTs. 
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APPENDIX A: PDF OF WORK OVER A TIME INTERVAL AND ENDPOINT 

POSITIONS 

Here, we will determine the Gaussian joint distributions P* of the work W T over a time interval 
of length r from time t to t + r and the positions of the Brownian particle at the beginning, Axt, 
and the end, Axt +T , where 

Ax t = x t -v*t. (Al) 



These are needed for the numerical sampling method in Sec. Mil as well as in the calculating of the 
Fourier transform of Pj 3 in Eq. (|36|) which is evaluated in App. [Bj We are interested both in the 
transient case, for which j = T and t = 0, and in the stationary state, for which j = S and t — > oo. 

For notational convenience, we introduce a seven dimensional vector a = (W T , Ax\, AX2). The 
PDF P* is characterized by the moments 



aj = JdaPj(a;T)a, (A2) 

which is a seven dimensional vector, and 

A = y'daP;(a;T)(a-a i )(a-a i ) t , (A3) 

Here the superscript dagger (f) denotes the transpose. It is clear from this definition that A is a 
real symmetric 7x7 matrix. Also, it has det A > 0. 
Once these moments are known, P* is given by: 

-i(a-ay)t-A-l.(a-&i) 

P* (a; r) = . ( A4) 

jK ' ^det(2^A) V ' 

Note that if det A = 0, then the PDF P* is a delta function (in one or more directions). 
To give the specific form for a and A, we first write 



/a: 



af I . (A5) 



(where is a scalar and a.^ and a.^ are three-vectors) and 



'An Al Al 

A = I A21 A 22 A^ 2 I , (A6) 

VA31 A32 A33, 



(where An is a scalar, A21 and A31 are three- vectors, A22, A23 and A33 are 3x3 matrices). 
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The specific form for a.j for the transient case, i.e., a T is obtained using Eqs. (JSJ), (fTHl) . (|3*U)) . 
PI) . dS3J, (US) and (dSJ, which yield 

4 1} = (VF T ) T = w(t - 1 + e" T ) (A7) 

4 2) = (Ax ) = (A8) 

4 3) = (Ax,) = (e~ T - l)v*. (A9) 

The sub-elements of A for the transient case are explicitly determined from Eqs. Q, Ijl0|) . (|1X|> . 
P|) : (EH), ifATj) . (IA21) and (|A5 |) - (|A9|) . giving 



An 


= ([W T - 


-(W T )] 2 ) T = 2w(T-l + e- T ) 


(A10) 


A 2 i 


= ([Ax 


-(Axo)][Wr-(Wr)]> 






= (e- T - 


l)v* 


(AH) 


A 3 i 


= ([Ax T 


-(Ax T )][W T -(W T )}) 






= (e- T - 


l)v* 


(A12) 


A 22 


= <[Ax 


- (Ax )][Ax -(Ax )] t > = l 


(A13) 


A33 


= ([Ax r 


- (Ax r )][Ax r - (Ax T )]t> = 1 


(A14) 


A32 


= ([Ax r 


-(Ax r )][Ax -(Ax )] t > = e - T l 


(A15) 



For the stationary case, the specific forms for the components of a s are found from Eqs. (jSJ), 
pi]). (|33l>. <EH>, fAljl . O and (EH): 

4 1} = {W T )s = WT (A16) 

4 2) = lim (Ax t ) = -v* (A17) 

af = lim (Ax t+r ) = -v* (A18) 

while the sub-elements of A are in that case, by Eqs. (EH1), (EH), W^ - and (JSEJ), 



An 


= lim ([W T - (W T )} 2 ) = ([W T - (W T )] 2 ) S 

t— >oo 


(A19) 


A21 


= lim ([Axt - <Ax t )][W T - (Wr)]) 


(A20) 


A31 


= lim ([Ax HT - (Ax t+T )][Ty r - (W T )}) 

t— >oo 


(A21) 


A22 


= lim ([Ax t - (Axt)][Axt - (Ax t )]t) 


(A22) 


A33 


= lim ([Axt +r - (Ax (+T )][Ax (+T - (Ax t+r )]t) 


(A23) 


A32 


= lim {[Ax t+T _(Ax t+T )][Ax t - (AxA]t) ; 


(A24) 



which turn out to be identical to those of the transient case in Eqs. lfA"T0jl - (fAT5} upon direct 
evaluation using Eqs. @-(HTJ), (EH), (|AH) and (|A16j) - (|A18|) . This is why we did not denote a j 
dependence of A. 
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APPENDIX B: FOURIER TRANSFORM OF tt? 

The Fourier transform of the PDF of heat will be calculated here, starting from Eq. IJ36JI . To 
calculate the Pj(q;r) from that equation, we define the quantities 

o ) (B1) 

'0 \ 

10, (B2) 
v -1/ 

so that one can write in the exponent in Eq. ()36l) as W T — \ {|Ax2| 2 — |Axi| 2 } = c • a + |a' • B • a. 
Using Eq. (|A4|) . one obtains then 

/„— |(a— aj)t •A _1 -(a— aj)+iiijat-B-a+i(jc-a 
rfa vmm • (B3) 

To evaluate this, the exponent is first rewritten as 

— -(a — a.j)^ ■ A^ 1 • (a — a.j) + 2*^ a ^ • B • a + iqc ■ a 

= - i(a - a^ • [A- 1 - iqB] ■ (a - a;.) + dj, (B4) 

where a'- = [I — iqA ■ B] _1 • (a,- + iq A • c) and 

dj = y (B • a., + c) f • (I - A • B)™ 1 • (a,- + iq A • c) + a,- • c . (B5) 

Here, I is the 7x7 identity matrix. Then substituting Eq. (|B4j) in Eq. ()B3|) and changing the 
integration variable to x = a — a!- yields 

Pj(q;r)= Le-^^A- 1 -^. ^ 
J 7 Vdet(2vrA) 

(B7) 



V^detO -igA- B)' 

where the identity det(A) det(B) = det(AB) has been used. To make Eq. (|B7|) into an explicit 
expression for Pj, the inverse of the matrix (I — iqA ■ B) is required in the expression for dj in 
Eq. ()B5|) . and in Eq. (|B7|) . its determinant. These are obtained as follows. Using Eqs. (|A6|) . 
(IXlOll - tlAToT) and (|B2|) . it follows that 

/l ig(e- T - l)v*t ig( e - T - l)v*t\ 
l-igAB= i<?e" r l . (B8) 

The determinant of this matrix is 

det(l - iqA • B) = [l + (1 — e" 2r )g 2 ] 3 . (B9) 
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For the inverse of Eq. I)B8|) , we get 



Q-iqA- B) 



ig(e- T -l)[l+ig(l-e~ T )]v*t tg(l- e - T )[l-tg(l-e- T )]v*t \ 
l+(l-e- 2T )q 2 l+(l-e-' 2T )q 2 



l+{l-e-' lT )q- 

l+(l-e-' 2T )q 2 A 
iqe~ T n 
l+(l-e-^)g 2 1 



(l-e- 2T )g 2 
l+(l-e- 2T ) 9 2 JL 
l+(l- e - 2T ) 9 2 1 



(BIO) 



/ 



We now have the material needed to calculate dj from Eq. I)B5|) explicitly. To calculate Eq. I)B5|) . we 
use Eqs. (|X5|) . (|Ail1|) . l|ATl^ . (jAT2"|) . fBT]) . (|H2|) and l|BT0|) . to find, after some rearrangements, 



l + (l-e 



-2t\„2 



iq 6 w{\ - e~ T f - iqe- T a.y> ■ a}° ; + ||a^|"(l + iq) - ±|a}° ; |"(l - iq) 



-r 5 (2) s (3) , l|-(2),2 



1 1^(3) ,2, 



ig(l - e~ r )[l + iq{l - e" T )]v* • af } + iq(l - e - T )[l - iq(l - e~ T )]v* • af } 



+ a- 1) + igw(r- l + e~ r )j 

Furthermore, from Eqs. (|A7|1 — (|A9|I and (|Bllj) . d T for the transient case follows as 

[l-e-][l + (i + 2g 2 )(l-e-)]\ 



(Bll) 



d T = wq(i — q) < r 



1 + (1 - e- 2r )g 2 



(B12) 



while using Eqs. (|A16|) - (|A18|) and (jBllj) . <i s for the stationary case, is, after some rewriting, found 
to be 



d s = wq(i - q) < r 



2<2 2 (1 -e 



-T\2 



(B13) 



l + g 2 (l-e- 2 -) 

Finally, by Eqs. ()B7|) and (|B9j) . these expressions for yield for the Fourier transforms explicitly 



PT(q;r) 
Ps(q;r) 



wq{i—q)* 

e 


L [l-e- T ][l + (l/2 + 2< ? 2 )(l-e-^)]l 


( l + (l-e- 2T )q 2 J 


[ 

wq(i-q) 

e 


1 + (1 - e" 2 -) 

2q 2 (l-e" T ) 2 


g2]3/2 


T l+(l-e-' 2T )q' 2 





[1 + (1 - e~ 2 -)g 2 ] 3 /2 



(B14) 
(B15) 
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